function [ bell_error,time] = bell_shape_error( position_x,position_y,velocity)

fs = 1000;
t = 0:1/1000:length(velocity)/1000;

% min jerk
out = min_jerk([position_x(1) position_y(1)], [position_x(end) position_y(end)], t);
bell_x = diff(out(:,1))*fs;
bell_y = diff(out(:,2))*fs;
bell_xy = sqrt((bell_x.^2+bell_y.^2));

% 10% peak bell
[max_v,position] = max(bell_xy);
onset_velxy = bell_xy(1:position);
%onset
IniNum = find(onset_velxy < 0.10 * max_v);
onset_bell = IniNum(end);
%offset
offset_velxy = bell_xy(position:length(bell_xy));
LastNum = find(offset_velxy < 0.10 * max_v);
offset_bell = LastNum(1) + position -1;

% 10% velocity
[max_v_1,position_1] = max(velocity);
onset_velxy_1 = velocity(1:position_1);
%offset
offset_velxy_1 = velocity(position_1:end);
LastNum_1 = find(offset_velxy_1 < 0.10 * max_v_1);
offset = LastNum_1(1) + position_1 -1;

valid_vel = velocity(1:offset);
time = length(valid_vel)/fs;
%resample
increment = length(valid_vel)/(offset_bell - onset_bell);
bell_shape = resample(bell_xy(onset_bell:offset_bell),fix(increment*fs),fs);

% RMSe = sqrt(sum((vel-bell).^2)/n)
[bell_rol,~] = size(bell_shape);
[velocity_rol,~] = size(valid_vel);

[min_rol,~] = min([bell_rol,velocity_rol]);
bell_error = sqrt(sum((bell_shape(1:min_rol-1)-valid_vel(1:min_rol-1)).^2)/(min_rol-1));
% figure(1)
% hold on,
% plot(bell_shape(1:min_rol-1),'color','r');
% plot(valid_vel(1:min_rol-1),'color','g');
% close;
